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We present a cluster Gutzwiller mean-field study for ground states and time-evolution dynamics in 
the Bose-Hubbard ladder (BHL), which can be realized by loading Bose atoms in double-well optical 
lattices. In our cluster mean-held approach, we treat each double-well unit of two lattice sites as 
a coherent whole for composing the cluster Gutzwiller ansatz, which may remain some residual 
correlations in each two-site unit. For a unbiased BHL, in addition to conventional superhuid phase 
and integer Mott insulator phases, we hnd that there are exotic fractional insulator phases if the 
inter-chain tunneling is much stronger than the intra-chain one. The fractional insulator phases can 
not be found by using a conventional mean-held treatment based upon the single-site Gutzwiller 
ansatz. For a biased BHL, we hnd there appear single-atom tunneling and interaction blockade if 
the system is dominated by the interplay between the on-site interaction and the inter-chain bias. In 
the many-body Landau-Zener process, in which the inter-chain bias is linearly swept from negative 
to positive or vice versa, our numerical results are qualitatively consistent with the experimental 
observation [Nat. Phys. 7 , 61 (2011)]. Our cluster bosonic Gutzwiller treatment is of promising 
perspectives in exploring exotic quantum phases and time-evolution dynamics of bosonic particles 
in superlattices. 

PACS numbers: 03.75.Lm, 05.70.Fh, 67.25.dj, 02.70.-c 


I. INTRODUCTION 

The unprecedented experimental techniques of ma¬ 
nipulating and detecting ultracold atoms in optical lat¬ 
tices [I|, 0 provide an ideal testing ground to investi¬ 
gate Bose-Hubbard (BH) models The remark¬ 

able cleanness and high tunability of ultracold atomic 
systems allow one to explore various many-body quan¬ 
tum phenomena in BH models dm. For an example, 
the experimental realization of the one-dimensional (ID) 
atomic Hubbard model [l^ provides new opportunities 
to exploring quantum statistical effects and strong cor¬ 
relation effects in low-dimensional quantum many-body 
systems [l^ . Quantum dynamics as well as quantum 
phase transition between superhuid (SF) phase and Mott 
insulator (MI) phase in BH models are of great interests 
and have been widely investigated [Hill- 

In recent, by loading ultracold Bose atoms into a 
double-well optical lattice potential, the Bose-Hubbard 
ladder (BHL) had been realized and the many- body 
Landau-Zener (LZ) dynamics has been explored (^ . 
Different from the single-particle LZ process, the break¬ 
down of adiabaticity in the inverse sweeping from the 
highest excited state had been observed in the many- 
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body LZ process of the BH ladder. This experiment has 
stimulated extensive investigation of both stationary and 
dynamic behaviors in the BHL via different theoretical 
methods, such as, full diagonalization method (2^ and 
time-dependent density-matrix renormalization group {t- 
DMRG) technique [24 - 1^ . However, the full diagonal¬ 
ization and t-DMRG methods should cost a huge number 
of computational resources. 

To simulate the BHL with less computational re¬ 
sources, the bosonic Gutzwiller method is an al¬ 

ternative option. Although the Gutzwiller method has 
common restrictions of the mean-field methods, it has 
provided versatile applications in qualitative calculations 
of both stationary states and time-evolution dynam¬ 
ics. In recent, cluster bosonic Gutzwiller methods [2^ 
Is^j have been developed by coupling multi-site clusters 
rather than single sites with the mean field. By em¬ 
ploying the cluster Gutzwiller method, some properties 
of many-body LZ phenomena in repulsive BHL [35j | and 
some quantum phases in attractive BHL have been 
explored. In Ref. [s^, the phase diagram and LZ dy¬ 
namics for fixed average number of particles per site have 
been shown. However, it does not discuss how the phase 
diagram and LZ dynamics depend on the chemical poten¬ 
tial. In the large-size multi-site Gutzwiller method [s^, 
the three-body constraint has been imposed to each lat¬ 
tice site. In a realistic experimental system, the num¬ 
ber of particles in each lattice site may break this con- 
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straint. In addition, the ground-state phase diagrams of 
BH systems have been obtained by the analytical mean- 
field approach , the cell strong-coupling perturbation 
techni que 1^ and the composite boson mean-field the¬ 
ory potklj etc. 

In this article, we present a cluster Gutzwiller mean- 
field study for the ground-state phase diagram and many- 
body LZ dynamics of a BHL. In our mean-field treat¬ 
ment, we regard each double-well unit of two lattice sites 
as a coherent whole for composing the cluster Gutzwiller 
ansatz, which remains some residual inter-site correla¬ 
tions in each double-well unit. For a unbiased BHL, in 
addition to superfiuid (SF) and integer Mott insulator 
(MI) phases which may be found by single-site Gutzwiller 
treatment, we find that there exist exotic fractional insu¬ 
lator phases if the inter-chain tunneling is much stronger 
than the intra-chain one. The exotic fractional insulator 
phase in BHL is similar to the rung-Mott phase in hard¬ 
core BH system [H, [s^ • We also obtain the phase dia¬ 
gram for the asymmetry BHL, which have not yet been 
reported. In further, by linearly sweeping the inter-chain 
bias from negative to positive or vice versa, we analyze 
many-body LZ dynamics in the system and confirm the 
existence of adiabaticity breakdown. 

This article is organized as follows. In Sec. im we give 
the physical model and discuss its realization. In Sec. lHIl 
we present the cluster Gutzwiller mean-field method and 
obtain the ground-state phase diagram for both sym¬ 
metric and asymmetric BHLs. In Sec. IlYl we study 
the many-body LZ dynamics and show the adiabaticity 
breakdown. At last, in Sec. El we briefly summarize and 
discuss our results. 


II. MODEL 

We consider an ensemble of Bose atoms confined within 
a double-well superlattice potential, 

V{x,z) = 14; sm^{27rx/Xa;i) -I- I4s sin^(27ra:/Axs) 

-I-14 sin^(27rz/A^), (1) 

where the first and second terms are generated by super¬ 
imposing two standing-wave lasers along the x-direction 
with wavelengths Xxi and Xxs ■ The two potential depths 
14s Enid 14; are determined by the laser intensities. To 
form the double-well lattice potential, the wavelengths 
are set to be X^i = 2Xxs- The last term describes a lat¬ 
tice potential along the z-direction with the wavelength 
A^ and the depth 14. During the experiment, the energy 
difference between the lattices in each double-well unit 
can be ramped up or down with time. The schematic 
diagram for the double-well lattices is shown in Fig. [T] 

If the barriers between neighboring double-well units 
along the x-direction is sufficiently high, the system can 
be described by several parallel BHLs with ignorable 
inter-ladder couplings. The Hamiltonian for a single BHL 



FIG. 1. Schematic diagram of the superlattice potential m 
projected onto the x-z plane. The time-dependent bias A(t) 
can be achieved by ramping up the lattices to obtain the tilted 
potential along the x-axis. Here, the parameters are set as 
A 2 = Xxi = 2Xxs and Vz = 1.5Ki = 1.5Ks. The poten¬ 
tial along the x-axis within a period of X^i is an asymmetric 
double-well potential and the potential along the z-axis is a 
standing-wave potential. The two lattice sites in each double¬ 
well potential are packed as a cluster, which are depicted by 
the gray rectangles. The clusters are decoupled by apply¬ 
ing the Gutzwiller mean-field treatment, in which the crosses 
stand for the decoupling between neighboring clusters and the 
gray dashed and solid lines respectively denote the intra- and 
inter-chain tunneling. 


reads as. 

Hit) = - J|| ^ - Jy ^ (bl^b.R + h.C.) 

{jk)<j j 

j 

— ^ ^ ' ffjo-; (2) 

where, (jk) indicates the summation comprising all near¬ 
est neighboring sites in the same chain and the index 
cr = iL,R) denotes the left or right chain. The symbols 
b^j^ (bja) creates (annihilates) a Bose atom on the j-th 
lattice of the cr-chain, and fij^ = b^j^bj„- stands for the 
atomic number. The parameters J|| and J± are the intra- 
and inter-chain nearest-neighbor hopping strengths, re¬ 
spectively. The on-site interaction U is determined by 
the s-wave scattering lengthes and the chemical poten¬ 
tial /i determines the particle filling. 

To characterize different regimes of the BHL, we intro¬ 
duce the ratio between the intra- and inter-chain hopping 
strengths /3 = J||/J_l. If /? <C I, the intra-chain tunnel¬ 
ing is ignorable and the ladder system can be regarded as 
isolated double-wells. However, if /3 ^ I, the intra-chain 
tunneling becomes dominant, the system can be treated 
as two decoupled single BH chains. The parameter A(t) 
represents the inter-chain energy bias. 
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III. GROUND-STATE PHASE DIAGRAM 

In this section, we show how to obtain the ground- 
state phase diagram via the cluster Gutzwiller mean-field 
treatment. In the first subsection, we describe the clus¬ 
ter Gutzwiller mean-field approach for the BHL. Then 
in the second subsection, we present the self-consistent 
procedure for determining ground states. In the last sub¬ 
section, we give the ground-state phase diagram. 


A. Cluster Gutzwiller mean-field approach 

The standard Gutzwiller method assumes the wave- 
function of the whole system as a product state of 
single-site wavefunctions. By implementing the standard 
Gutzwiller procedure, the BH model is decoupled as sin¬ 
gle sites which couple with surround sites via their aver¬ 
age mean fields. In further, attribute to the equivalence 
of all lattice sites in the model, one can replace the mean 
fields of surround sites with the mean field of the site 
itself and so that the mean-field version for the original 
Hamiltonian can be written as a sum of single-site terms. 

In the following, the cluster Gutzwiller mean-field ap¬ 
proach is an extension of the single-site Gutzwiller mean- 
field approach. As shown in Fig. [U the decoupling holds 
for each double-well cluster which includes one lattice site 
in the left chain and one lattice site in right chain. There¬ 
fore all clusters are equivalent and the state for the whole 
system is written as a product state of the single-cluster 
states which remains the correlations between lattice sites 
in the same cluster. Unlike the single-site Gutzwiller ap¬ 
proach, in which all tunneling terms are decoupled, the 
cluster Gutzwiller approach keeps the intra-cluster tun¬ 
neling terms and only decouple the inter-cluster tunnel¬ 
ing terms. By using the mean-field treatment, the inter¬ 
cluster tunneling terms are decoupled as 

blh, « blik.) + (bl)k. - {bl){ka) 

= b]a^ka + ^P%bkc7 - (3) 

where Lpka = {bka) and the high-order fluctuations 
= (kjcr-{kj<y)){ka-{ka)) are neglected. There¬ 
fore the original Hamiltonian ([2]) is decoupled as 

(4) 

3 

with the single-cluster mean-field Hamiltonian 
^MF ^ ^ {^^kab]^ + ^labja - Re [p^a^fka] ) 

— J_L (kjj^bjJi -f h.c.^ -)- — flja {pja — i) 

(7 

(5) 


Making use of the Gutzwiller ansatz, the state for the 
whole system can be expressed as a product state of 
single-cluster states, 

ivi,GA)^ni'k)^., ( 6 ) 

3 

where the state for the j-th cluster \'^)j can be expanded 
as 

\'^h=J2 E fNl\N,m)^. (7) 

m——N 


with I TV, m)^ denoting the state basis for the j-th cluster. 
Here, TV = TV^ -|- Nr, m = Nr — Nr, Nr {Nr) stands 
for the number of particles in the left (right) chain, the 
probability amplitudes are complex numbers, and 

TVniax is the truncation of the maximum particle number. 

Obviously, it is easy to find that the eigenequation 
^MF jxjfGA^ _ |v[/GA^ Jqj. whole System is equiv¬ 

alent to the single-cluster eigenequation 

= (8) 


with E = Ej ■ By substituting the single-cluster 
state © into the single-cluster eigenequation ([5|) , we ob¬ 
tain 


^3JN,m ~ 


J|| 


- ^(l)jRy/N + 


J\\ 




J\\ 


cU) 






AT+l ,m— 1 

U) 


•/||Re [p*jr(I> 3L + P*3R(t>iR] In]. 
- iyVTV-fm\/TV-TO-b 

J± 

2 


— —E^N + m + 2\/TV — 


N,m-\-2 


^ (TV^ -|- — 2TV) -I- —m — fj,N 


fU) 

JN,m- 

(9) 


Here, the order parameters are quantum expectation 
values of bosonic annihilation operators, i.e. (pjR = 
{'i>^^\bjR\'^^^) and pjR = After some 

mathematical calculation, we have 


N + m + 2 (j) 

PjL — 2^ Y 2 TAr,mTAr+l,m-|-H 
N,m 



N — m + 2 
2 


fU)* fU) 

N^m3 N-\-l^m— 


1- 


( 10 ) 

( 11 ) 


For convenience, we define (pjR = + p(^j_i')R and 

4>3R = '^{3 + ^)R + ^(3-i)R- 
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B. Self-consistent procedure for determining 
ground states 

As the single-cluster Hamiltonian ([5]) depends on the 
mean fields, one has to implement self-consistent pro¬ 
cedure for determining the mean fields and the ground 
states. Given the parameters U, Jy, Jj_, /i and A, one can 
obtain the ground state from the single-cluster eigenequa- 
tion (|8l) and the self-consistent relations (pa = (ba)- In 
Fig. 2, we show the key steps of the self-consistent pro¬ 
cedure for determining ground states. 

(i) Initialize ipcr = 0, ‘p'a = ‘fa nnd set the trial ground 
state energy £’^ 5 " an arbitrary value. 

(ii) Substitute ipcr into the single-cluster Hamiltonian, 
and diagonalize the Hamiltonian to obtain its ground 
state with eigen-energy Ecs- 

(iii) If Egs < replace and with Egs 

and ipa-, respectively. Otherwise, let ipa- = fa + ^fa and 
implement step (ii) again. 

(iv) Repeat Steps (ii) and (iii) until fa- > VNmax ■ 

(v) Set fa = f'a and calculate the ground state |GSi) 
from the single-cluster eigenequation. 

(vi) Calculate the order parameters f^a lor the ground 
state |GSi). 

(vii) Compare f]a and if \f]j — f^a\ < e (where e 
is pre-given tolerance), then output |GSi) as the ground 
state. Otherwise, set fa = f]} and return to step (v). 

Through the procedure from step (i) to step (iv), one 
can numerically minimize the system energy with re¬ 
spect to the order parameters in the interval of fa € 
[0, VA^max]. Usually, the steps (v-vii) are the so-called 
self-consistent procedure. 


For parameters and tolerance: e 


initial = 0, and set cp'^ = (p^ 


diagonal {(p^) 
find the ground energy: £, 



diagonal 
find the ground state: | GS,) 
^ ^:=(GS,|i|GS,) 



FIG. 2. The numerical self-consistent procedure for solving 
the single-cluster eigenequation. 


C. Phase diagram 

In this subsection, we show the ground-state phase di¬ 
agram for the BHL. Our cluster mean-field approach can 
be applied to both symmetric and asymmetric systems. 
Below, we first consider the symmetric system with no 
inter-chain bias (i.e. A(t) = 0), then consider the asym¬ 
metric cases with nonzero inter-chain bias A. 

In Fig. [21 we show the ground-state phase diagram for 
symmetric BHL with A(t) =0 and different ratios /3 = 
J||/J_L. Due to the absence of asymmetry, the two chains 
are completely equivalent and the order parameters of 
both chains are always equal fjL = fjR = fj- Therefore, 
it is enough to give the phase diagram via analyzing the 
order parameter for one of the two chains. The ground 
states sensitively depend on the chemical potential /r, the 
on-site interaction U, the intra-chain hopping Jy and the 
inter-chain hopping J_l. 

Usually, determined by the superfluid order parame¬ 
ter, the BH systems have two typical phases: (i) the 
superfluid (SF) phase of nonzero order parameter and 
(ii) the Mott insulator (MI) phase of zero order param¬ 
eter and integer filling number. For our atomic BHL 
system, when the inter-cluster hopping Jy is sufficiently 


strong, the atoms can move freely between neighboring 
double-well clusters and there appears a SF along the 
chain direction. In contrast, when the on-site interaction 
U becomes sufficiently strong, the atoms are localized in 
each cluster and there is no SF along the chain direction. 
The chemical potential p controls the filling number, i.e. 
the average atomic number per site. 

Under the condition of /3 = J\\/J± ^ 1, i.e. the intra¬ 
chain tunneling is much stronger than the inter-chain 
tunneling, the BHL can be regarded as two decoupled 
chains and the corresponding phase diagram is almost as 
same as the one for a single BH chain. In Fig. [2] (c), 
we show the phase diagram for /3 = 10. At the side of 
strong intra-chain tunneling, J\\/U —>■ -l-oo, the ground 
states are SF phases of nonzero order parameter. At the 
side of strong interaction, J\\/U —>■ 0, there appear sev¬ 
eral integer MI lobes which has integer filling numbers per 
lattice site and zero order parameter. The blue region in 
the bottom corresponds to the vacuum state with no any 
atoms. The biggest lobe corresponds to the MI phase 
with definitely one atom (n = 1 ) in each site and the 
smaller one stands for the MI phase of n = 2. This phase 
diagram reminds us the one for the one-dimensional BH 
model (il ]. 
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FIG. 3. (Color online) The ground-state phase diagram for the symmetric Bose-Hubbard ladder (A = 0) with the on-site 
interaction [7 = 1 and different values oi P = J||/Jx: (a) /3 = 0.1, (b) P = 1, and (c) P = 10. In our calculation, we set 
the truncation of maximum particle number Amax = 6 whose validity has been numerically verified. The blue areas are the 
insulator phases with zero order parameter ipj = 0 with n denoting the filling number (the average atomic number per lattice 
site). The Mott insulator (MI) lobes have integer filling numbers, while the loophole insulator (LI) phases have half-integer 
filling numbers. The regions outside the blue areas are the superfluid (SF) phases with nonzero order parameters ipj 7 ^ 0. 


The areas of MI lobes shrink if the ratio /3 decreases, 
see Fig. [3] (a-c). Qualitatively, the shrinking of MI lobes 
can be understood by the intra-chain tunneling assisted 
by the inter-chain tunneling. When the ratio /3 becomes 
very small, the inter-chain hopping J± are much stronger 
than the intra-chain hopping J||, the areas of MI lobes 
shrink dramatically, and, interestingly, several loophole 
insulator (LI) phases of zero order parameters appear 
between the conventional MI lobes, see Fig. |3](a). 

To distinguish the LI and MI phases, we calculate the 
filling numbers (the average atomic numbers per site) 
and find that the LI phases have half-integer filling num¬ 
bers and while the MI phases have integer filling num¬ 
bers. The half-integer filling numbers mean that the to¬ 
tal atomic numbers per cluster are odd integer numbers 
and the residual atom in each double-well cluster can 
freely move between the two wells of each cluster. In 
further, we calculate the intra-cluster first-order correla¬ 


tion Cor|(^^ = 



and find that the LI phases have 


nonzero Cor|j^^. 

The appearance of the LI phases is a direct result of 
U ^ J|| and J± J||. As U J||, the tunneling along 
the chain direction is suppressed and the order param¬ 
eter vanishes. However, the atoms in each double-well 
cluster may still freely move between the two wells and 
so that the total atomic numbers per cluster may be 
odd integer numbers. The insulator phases of fractional 
filling numbers have also been found in one-dimensional 
superlattice BH models [H, via mean-field method, 
quantum Monte Carlo simulation and numerical density 
matrix renormalization group simulation. Different from 
the one-dimensional superlattice BH chains [l^, , our 


ladder system includes two coupled one-dimensional BH 
chains and the coupling between different clusters are 
more complex. 

Now, we discuss the ground-state phase diagrams for 
BHL with nonzero bias A. Due to nonzero bias A, the 
two chains are no longer equivalent and so that the order 
parameters ipjL and ipjn for the left and right chains may 
have different values. In Fig. |4l we show the two order 
parameters {<PjL, ^jr) (the second and third rows) and 
the intra-cluster first-order correlation Cor||^^ (the first 
row). 

In our numerical simulation, by employing the cluster 
mean-field method presented in Subsection A, we con¬ 
sistently obtain the ground-state phase diagram in the 
(A/f7, ^/f7)-plane for the intra-chain hopping strength 
J|| = 0.01, the on-site interaction U = 1, and different 
values of the hopping ratio P = J\\/ J± = (0.1,1,10). Our 
results show that, for given Jy, /3 and the order param¬ 
eter for the left chain with bias A equals to the order 
parameter for the right chain ipjfi with bias —A. There¬ 
fore, the ground-state phase diagrams of (pjL and ipju 
are symmetric with each other about the axis A = 0 for 
given J|| and p. Moreover, the intra-cluster correlation 

Cor||^^ is also symmetric about the axis A = 0. 

If the inter-chain tunneling is very weak, that is the 
hopping ratio P ^ 1, the BHL can almost be treated as 
two independent chains. In the third column of Fig. 4, 
we show the inter-chain coherence Cor||^^ and the order 
parameters {ipjL, ‘fijn) for P = 10. In phase diagrams 
of {^jLi VjR)^ there appear several parallel and equal¬ 
spaced SF stripes with nonzero order parameters ipjL or 
iPjH, see Fig. 4 (f) and (i). The SF stripes of ipjL become 
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FIG. 4. (Color online) The ground-state phase diagrams for the biased Bose-Hubbard ladder with different values of /3 = J\\/ J_l 
versus the bias A. In our calculation, we fix the on-site interaction U = 1, the intra-chain hopping Jy = 0.01, and set the 
truncation of maximum particle number Amax = 6. The first row shows the intra-cluster first-order correlation Cor^^, and the 
second and third rows respectively correspond to the order parameters <Pl,r for the left and right chains. The hopping ratios 
are chosen as /I = 0.1 (the first column), /3 = 1 (the second column) and /3 = 10 (the third column), respectively. 


narrower and its corresponding values become smaller as 
the bias A increases from the negative to the positive 
side, see Fig. 4 (f). While for the order parameter ipjR^ 
its SF stripes change oppositely as they are symmetric 
with the ones of the order parameter about the axis 
A = 0, see Fig. 4 (i). The blue regions correspond to MI 
phases with zero order parameters ipjL or <fjR. For the 
inter-chain coherence Cor||^\ nonzero values only appear 
in the vicinity surrounding some specific points, which 
form an inverted triangle lattice structure, see Fig. 4 (c). 

The interplay between the inter-chain bias and the on¬ 
site interaction under weak hopping leads to the exotic 
single-atom tunneling [the bright spots in Fig. 4 (c)] and 
the interaction blockade [the blue regions in Fig. 4 (c)]. 


Under strong on-site interactions, that is U ^ Jy and 
U ^ J±, the system behavior can be understood by a 
perturbation picture. As C/ » Jy and U J_\_, the 
hopping terms can be regarded as perturbations and the 
dominant part of the BHL Hamiltonian reads as 

Hd = — nja {pjrr ~ l) ^ ~ ^1^) 

fo- j 

( 12 ) 

fo- 

From the Fock states for a cluster, a single Fock state 
{til, Ur) corresponds to a MI phase of zero order param¬ 
eters, the quasi-degeneracy of \nL,nR) and jna -I- IjUr) 
will result a nonzero order parameter ipjR , and the quasi- 
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degeneracy of \ni,,nf^ and + 1 ) will induce a 

nonzero order parameter ifiR- 

In the MI regions of (pjL (i-e. ipjL = 0), the clus¬ 
ter state is a single Fock state \nL,nR). The quasi¬ 
degeneracy of I riL, ) and | nz, -I-1, n/j) under weak intra¬ 
chain hopping, which means that one atom can freely 
move in the left chain, will result a nonzero order pa¬ 
rameter ifR. From the energy quasi-degeneracy relation 
E{nL,nii) = E(nL + l^nu), we have the system param¬ 
eters obeying 

y-hC/ni,(ni =0,1,2,...). (13) 

Obviously, the above relation (US between the chem¬ 
ical potential /i and the bias A well agree with the SF 
stripes shown Fig. 4 (f). Similarly, from the energy quasi¬ 
degeneracy relation E{nL, tir) = E^ur, ur -I- 1 ), we have 

fi = + UnR,{nR = 0,1,2,...), (14) 

for the SF stripes of nonzero order parameter (pjR shown 
in Fig. 4 (i). 

In addition to the single-atom tunneling along a spe¬ 
cific chain, there exists single-atom tunneling between 
two chains, which corresponds to a nonzero intra-cluster 
first-order correlation Cor||^\ The inter-chain single¬ 
atom tunneling is caused by the quasi-degeneracy of 
\nL,nR) and \nL + ^,nR—\). Therefore, from the energy 
quasi-degeneracy £'(ni, uzj) = E(nL + l,nR — l), we have 

^ = U{nR-nL),{nL,R = 0,1,2,...) (15) 

As nonzero Cor^^ appears in the region of both (pjR ^ 0 
and (fjR 7 ^ 0 , the corresponding chemical potential /r is 
given as. 


^J-= l^{nL + nR), (16) 

with riL.fi = (0,1,2,...). This means that nonzero Cor^^ 
appears in the vicinity surrounding {A*/U,fi*/U) = 
{ur - Ur, ^{ur + Ur)) with ur^r = (0,1, 2,...), See the 
bright spots in Fig. 4 (b) and (c). For a given chemical 
potential fi = ^{ur + nR)U, a sequence of single-atom 
tunneling and interaction blockade takes place when the 
bias A increases from negative infinity to positive infin¬ 
ity. The single-atom tunneling and interaction blockade 
in the BHL with fixed chemical potential is reminiscent 
of that of the quantized Bose-Josephson junction with 
strong interaction To the best of our knowledge, 

the single-atom tunneling and interaction blockade in the 
BHL have never been reported before. 

In the second column of Fig. 4, we show the phase 
diagrams for the case of /? = 1. The parallel SF stripes of 
(fijR or ipjR still appear but blur at the quasi-degenerate 
regions in the vicinity of {A*/U, /U). Different from 

the case of large /3, (pjR ^ 0 and i^^r ^ 0 may coexist 
in some specific regions. Correspondingly, due to the 


increase of Jr, the area of nonzero Cor|j^^ surrounding 
{A* jU, fi* jU) extend. 

In the first column of Fig. 4, we show the phase di¬ 
agrams for the case of /3 = 0.1. The strong intra¬ 
chain tunneling makes the occurrence of inter-chain intra¬ 
cluster single-atom tunneling more easy, the regions of 
Cor||^^ 7 ^ 0 extend and merge into an entire area. Cor¬ 
respondingly, due to the strong inter-chain hopping, the 
properties of ipjR and ipjR change dramatically. The par¬ 
allel SF stripes are tailored and several avoided crossings 
appear in the vicinity of {A*/U, fi*/U). The avoided 

crossings, which have Cor|^^ 7 ^ 0 and zero order parame¬ 
ters {^jR = ipjR = 0), correspond to the LI phases shown 
in Fig. 3 (a). This means that atoms may move freely 
between the two chains although there is no superfluid 
along the chain direction. 


IV. LANDAU-ZENER DYNAMICS 


In this section, we analyze the many-body LZ dynam¬ 
ics in the BHL. In the many-body LZ process, the inter¬ 
chain energy bias A{t) is linearly swept from negative to 
positive or vice versa. The linear sweep of bias is de¬ 
scribed by A(f) = Aq + at with Aq being the initial bias 
and a denoting the sweeping rate. In the first subsec¬ 
tion, we show how to apply the Gutzwiller mean-field to 
the time-evolution problem of our BHL system. In the 
second subsection, we present the population dynamics 
in the ground-state sweep and the inverse sweep, respec¬ 
tively. In the ground-state sweep, the initial state is the 
ground state, while in the inverse sweep, the initial state 
is the highest excited state. 



FIG. 5. The numerical simulation procedure for the time- 
evolution of Bose-Hubbard ladder system via our cluster 
Gutzwiller mean-field method. 















A. Time-evolution problem 

The time-evolution obeys the Schrodinger equation 

(17) 

where H {t) is the original time-dependent Hamilto¬ 
nian By applying the dynamical Gutzwiller mean- 
field method, the time-evolution is described by the dy¬ 
namical Gutzwiller equations 

= (18) 

where is the time-dependent mean-field Hamil¬ 

tonian and |'I''^^(t)) = r[^. l'I'(t))^- denotes the time- 
dependent Gutzwiller ansatz. Here, the single-cluster 
state reads as 

-A/^max N 

E E fN]mit)\N.m)r ( 19 ) 

m——N I 


Similar to determining the ground states, the time- 
dependent mean-field Hamiltonian can be decou¬ 

pled as a sum of single-cluster Hamiltonians. Therefore, 
the dynamical Gutzwiller equations (fT51) can be simpli¬ 
fied to the single-cluster equations 

= ( 20 ) 
with the time-dependent single-cluster Hamiltonian 
ijMF(t) = - J|| (^V’kabl^ + ^labja “ Re [<P*^^ka] ) 

(T,k=j±l 

— J± {b^jlPiR + h-C.^ -f — {p’jcr — i) 

a 

2 (^it? ~ RjL) ~ M ^ ( Rja, (21) 

a 

in which the time-dependent order parameters are given 
as = {'^j{t)\bja\'ij{t)). Substituting Eq. (fT^ and 

Eq. (I^Tl) into Eq. (I20p . one can obtain the following dif¬ 
ferential equations for the expansion coefficients 




J|| 


- -^(l)jL{t)VN + m/7V-l.m-l(t) - -^(j)jR{t)VN - mfN-l,m+l{t) 

- -^4>jL{t)*VN + m + 2fN+l,m+l{t) - -^(l)jR{t)*VN -m + 2fN+l,Tn-l(,t) 

- ^\/N + my/N -m + 2fN,m-2{t) - N + m + 2\/N - m/Ar,m+ 2 (t) 


U 


{N‘ 


2 I 2 

' m 


-27V) 


A(t) 


m- fiN + Ji|Re [(pjL{t)*(j)jL{t) + <fjR{t)*(j)jR{t)] 


fN,m{t), (22) 


with the time-dependent order parameters 

= 7’i-i-i,o'(l) T i,o-(t); (23) 

^jL{t) = Y y ^ ^ + ^ /ffm(0/^|l,m+l(l): (24) 

iV,m ^ 

‘fjRi.t) = E v/ ^ ~ ^ ^ /ffm(l)/^|l,m-l(0- (25) 

N,m 

By using the fourth-order Ronge-Kutta method, we sim¬ 
ulate the dynamics obeying Eq. (1221) . The flow chart for 
the numerical procedure is shown in Fig. [SJ Given the 
parameters J||, J_l, U, the initial bias Aq, the sweeping 
rate a, and the initial state, the time-dependent order pa¬ 
rameters should be estimated by the instantaneous states 
step by step. That is, for a specific time step, based upon 
the current state and the current order parameters, we 
need estimate not only the time-dependent state but also 
the time-dependent order parameters for the next time 
step. 


B. Population dynamics 

We consider two typical sweep processes: the ground 
state sweep and the inverse sweep. In the ground-state 
sweep, the system is prepared in the ground state of all 
particles in the lower chain, and the initial bias between 
left and right chains is set as Ag = —50 and then the 
bias A(t) is linearly swept from Aq to — Ag with the 
sweep rate a = —2Ag/r > 0. In the inverse sweep, 
the system is prepared in the highest excited state of all 
particles in the higher chain, and the initial bias between 
left and right chains is set as Ag = 50 and then the 
bias A(t) is linearly swept from Ag to — Ag with the 
sweep rate a = —2 Aq/T < 0. Here T is the total sweep 
time. For convenience, we assume that the initial state 
for both two sweep processes is the state of all atoms in 
the left chain. To show the many-body LZ dynamics, we 
calculate the transfer fraction nR{t) = NR{t)/N, which 
is the fraction of the particles in the right chain at a given 
time t. Obviously, the bias A(t) vanishes at time t = T /2, 
which corresponds to an instantaneous symmetric BHL. 
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Inverse Sweep 
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FIG. 6. (Color online) The many-body Landau-Zener dynamics in the Bose-Hubbard ladder for different sweeping rates |a|. 
Here, nR{t) stands for the transfer fraction, the cutoff of the maximum particle number is fixed as A^max=6. The other 
parameters are chosen as: |Ao| = 50, jj. = 3, U = 0.5, Jy = 0.25, J± — 1. Red-dashed lines and blue-solid lines correspond to 
the ground-state sweep and the inverse sweep, respectively, (a) For a large sweep rate, lal = 25, the transfer fractions n_R(T) 
for both the ground-state sweep and the inverse sweep are much smaller than 1. The difference between the the two sweeps is 
small, (b) For an intermediate sweep rate, |a| = 5, the transfer fractions nniT) for the ground-state sweep increases to almost 
1, and nii{T) for the inverse sweep is still far below 1. (c) For a small sweep rate, |a| = 1, the transfer fractions nniT) for 
the ground-state sweep reaches 1, while nR{T) for the inverse sweep is still below 1. This means that the ground-state sweep 
evolves adiabatically, while the inverse sweep evolves non-adiabatically. 


If there are no intra-chain hopping and no on-site in¬ 
teraction, i.e. J|| = 0 and U = 0, the physical picture for 
the many-body LZ dynamics is as same as the one for the 
conventional two-level LZ problem. This means, the final 
transfer efficiency is given by the conventional LZ formula 
ni{(-|-oo) = 1 — exp(—27rJ|/?i|a|) and there is no sig¬ 
nificant difference between the ground-state and inverse 
sweeps. However, taking into account the on-site inter¬ 
action and the intra-chain hopping, the many-body LZ 
dynamics becomes very different from the conventional 
two-level LZ problem. Below, we analyze the many- 
body LZ dynamics for the on-site interaction U = 0.5, 
the inter-chain hopping J_l = 1, the intra-chain hopping 
J|| = 0.25, and different sweep rates a. 

Independent on the sweep rate a, significant popula¬ 
tion transfers from the left chain to the right chain appear 
around the time t = T /2. This significant population 
transfer between the two chains is caused by the avoided 
level crossing in the vicinity of the bias A(t) = 0. How¬ 
ever, the transfer fraction nii{t) = NR{t)/N sensitively 
depends on the sweep rates, the physical parameters and 
the initial states. In particular, for slow sweep rates, 
there appear significant difference of the final transfer 
fraction for the ground-state sweep and the inverse sweep, 
see Fig. [51 

For a large sweep rate, |a| = 25, both the ground- 
state sweep (a = -1-25) and the inverse sweep {a = —25) 
are non-adiabatic, see Fig. [5] (a). The dynamics of the 
ground-state sweep and the inverse sweep is very simi¬ 
lar. The transfer fraction unit) rapidly increase around 
A(t) = 0 and then keep oscillates around a specific value. 
The final transfer fraction nR{T) is much below 1 because 


of the non-adiabatic evolution under fast sweeps. 

For an intermediate sweep rate, |a| = 5, the non- 
adiabatic excitation in the ground-state sweep is not very 
significant, while the non-adiabatic excitation in the in¬ 
verse sweep is very significant, see Fig. [n](b). After the 
system goes through the avoided level crossing region 
around A(t) = 0, the transfer fraction for the ground- 
state sweep (a = 4-5) is very close to 1 and its oscilla¬ 
tion amplitude is very small. While in the inverse sweep 
(a = —5), the final transfer fraction is much below 1 and 
the corresponding oscillation amplitude is much larger 
than the one for the ground-state sweep. 

For a small sweep rate, |a| = 1, the ground-state sweep 
undergoes adiabatic evolution and but the inverse sweep 
still show significant non-adiabatic excitations, see Fig. [5] 
(c). In the ground-state sweep (a = 4-1), there is no sig¬ 
nificant oscillations in the transfer fraction and the final 
transfer fraction is almost the perfect limit nji[T) = 1, 
which means that all particles in the left chain can be 
completely transferred into the right chain. However, in 
the inverse sweep (a = —1), the final transfer fraction is 
still much below 1 and the oscillation amplitude is still 
very significant, which indicates that there still exist sig¬ 
nificant non-adiabatic excitations. 

The adiabaticity breakdown in the inverse sweep quali¬ 
tatively explains the recent experimental observation [^, 
[^ . The observed adiabaticity breakdown, which can 
not be found in the conventional two-level LZ problem, 
is a result of the inter-particle interaction. Due to the 
inter- par ticle interaction, swallow-tail-shaped loop struc¬ 
tures |45l . l46l| , which correspond to the macroscopic quan¬ 
tum self-trapping in mean-field models (ttI - I^ . may ap- 
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pear in the energy spectrum for our BHL system. Unlike 
the conventional two-level LZ problem, whose energy- 
level structures for the ground state and the highest- 
excited state are similar, the energy-level structures for 
the ground state and the highest-excited state of our 
BHL system are very different. Because of their differ¬ 
ent energy-level structures, the ground-state sweep and 
the inverse sweep show different adiabatic/non-adiabatic 
dynamics. 

V. CONCLUSION AND DISCUSSION 

In summary, we present a cluster Gutzwiller mean-field 
approach to explore the static and dynamical behavior of 
the BHL, which can be experimentally realized by loading 
Bose atoms into a double-well optical superlattice poten¬ 
tial. In our mean-field treatment, the wavefunction of the 
whole system is assumed in form of the Gutzwiller ansatz, 
the two sites in each double-well unit are packed as a 
cluster and the inter-cluster hopping is decoupled by us¬ 
ing the conventional mean-field approximation. Through 
implementing the numerical self-consistent procedure, for 
both unbiased and biased BHLs, we obtain the ground 
states and give the phase diagram by calculating the or¬ 
der parameters. 

For an unbiased BHL, if the intra-chain hopping is 
much stronger than the inter-chain hopping (i.e. (3 « 
1 ), there appear several exotic loophole-shaped insula¬ 
tor regions of the half-integer filling numbers, which lie 
between the conventional MI lobes of integer filling num¬ 
bers. As (3 increases, the loophole-shaped insulator re¬ 
gions gradually shrink and disappear. Differently, if the 
inter-chain hopping is much stronger than the intra-chain 
hopping (i.e. (3 » 1), the unbiased BHL system can be 
regarded as two single BH chains and the corresponding 
phase diagram is almost as same as the one for a single 
BH chain. 

For a biased BHL, single-atom tunneling and inter¬ 
action blockade appear if the hopping terms are weak 
enough to be treated as perturbations. We present an an¬ 
alytical interpretation for the single-atom tunneling and 
interaction blockade based upon the quasi-degeneracy of 
different Fock states for the considered cluster. If the 
inter-chain hopping is much stronger than the intra-chain 
one, there appear exotic LI phases with no superfluids 


along the chain direction but nonzero inter-chain coher¬ 
ence. 

In further, we analyze the many-body LZ process of 
the BHL, in which the inter-chain bias is linearly swept 
from positive to negative or vice versa. We consider two 
different sweeps: the ground-state sweep and the inverse 
sweep. In the ground-state sweep, the initial state is 
the ground state and the final transfer fraction can reach 
1 if the sweep rate is small enough. While in the in¬ 
verse sweep, whose initial state is the highest excited 
state, there still exist significant non-adiabatic excita¬ 
tions when the corresponding ground-state sweep obeys 
adiabatic evolution. The breakdown of adiabaticity in 
the inverse sweep, which are well consistent with the 
recent experimental observations [l^, [2^ , is a result of 
the swallow-tail-sha ped loop structures induced by inter¬ 
particle interaction |45l . . 

In recent, for ultracold atoms in optical lattices, arti¬ 
ficial gauge fields have been realized by lattice shaking 
technique or laser-induced tunneling The ar¬ 

tificial gauge fields, which allow one to generate spin- 
orbit couplings and effective magnetic fields, opens a 
new path to explore quantum Hall effect and topolog¬ 
ical phases of matters. Our cluster Gutzwiller mean- 
field approach can also be extended to investigate the 
bosonic ladders in the presence of an artificial magnetic 
field 1 ^. IstI - I^ . such as the observation of chiral cur¬ 
rents |57l |. the measurement of Chern number in Hofs- 
tadter bands [H, , and the two-leg Bose-Hubbard lad¬ 
der under a magnetic flux . In addition, our clus¬ 

ter Gutzwiller mean-field approach may also use to ex¬ 
plore the non-equilibrium dynamics of two coupled one¬ 
dimensional Luttinger liquids and the dynamical in¬ 
stability of interacting bosons in disordered lattices [i^ . 
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